Power-law friction in closely-packed granular materials 
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In order to understand the nature of friction in dense granular materials, a discrete element 
simulation on granular layers subjected to isobaric plain shear is performed. It is found that the 
friction coefficient increases as the power of the shear rate, the exponent of which does not depend 
on the material constants. Using a nondimensional parameter that is known as the inertial number, 
the power law can be cast in a generalized form so that the friction coefficients at different confining 
pressures collapse on the same curve. We show that the volume fraction also obeys a power law. 



Friction is one of the oldest problems in science because 
it dominates various phenomena in our daily life. In par- 
ticular, dynamics of granular flow, which is ubiquitous 
in earth sciences and engineering, is governed by a law 
that describes behavior of the friction coefficient (ratio of 
the shear stress to the normal stress) . Such examples are 
avalanche, landslide, debris flow, silo flow, etc. In addi- 
tion, the nature of friction on faults, which plays a key 
role in earthquake mechanics [l|, , is also attributed to 
that of granular rock because fault zone consists of lay- 
ers of fine rock particles that are ground-up by the fault 
motion of the past. To find a suitable law of friction in 
granular materials under a specific condition is thus an 
essential problem. 

Although the frictional properties of granular materi- 
als are so important, our understanding is still limited. 
In the context of earthquake mechanics, slip velocity (or 
shear rate) dependence of friction coefficient, which is 
equivalent to rheology under constant pressure condition, 
is a matter of focus [2[ . In experiments on thin granular 
layers that are sheared at relatively low sliding velocities 
ranging from nm/s to mm/s, the behavior of the friction 
coefficient can be described by a phenomenological law 
in which friction coefficient logarithmically depends on 
the sliding velocity. This is known as the rate and state 
dependent friction (RSF) law ||. Note that the RSF law 
also applies to friction at interfaces between two solids, 
as well as that in granular layers. Although the RSF law 
applies well to lower speed (creep-like) friction, it is vio- 
lated in high speed friction. For example, several exper- 
iments indicated nonlogarithmic increase of the friction 
coefficient in granular layers at higher sliding velocities 
0) IE Hi- The same tendency was also observed in ex- 
periments on friction between two sheets of paper @, [1] . 
However, at this point, we do not know any friction law 
that is valid at such higher velocities. 

Several recent attempts to understand the nature of 
friction in granular media under high shear rates are note- 
worthy here. Jop and the coworkers presented a simple 
friction law that describes flow on inclined planes Q, 
based on massive simulations and experiments [Tfl |. Al- 
though their friction law seems feasible, it involves rather 
dilute flow and its applicability to denser and slower flow 
(e.g. quasistatic flow) is not clear. Da Cruz et al. 
performed an extensive simulation that focused on dense 
and slow regime and found a friction law that does not 



contradict that of Jop et al. However, because da Cruz 
et al. involved a two dimensional system, effect of the 
dimensionality may be questioned. In particular, in qua- 
sistatic regime where the nature of interparticle contacts 
plays an esential role in rheology, the effect of dimension- 
ality should be seriously investigated. 

In this Rapid Communication, we perform a three di- 
mensional simulation in order to understand the nature 
of friction in a slowly-sheared dense granular material. 
Our particular interest is a dense granular matter un- 
der high confining pressure (e.g. tens of MPa) which is 
roughly corresponds to a typical configuration of faults 
at seismic slip. Note that the RSF law is violated in 
such a situation. A new law is reported in which friction 
coefficient increases as the power of shear rate. 

In the following we describe the computational model 
of granular layers. The individual constituents are as- 
sumed to be spheres, and their diameters range uni- 
formly from 0.7d to l.Od. The interaction force follows 
the discrete element method (DEM) [!§]■ Consider a 
grain i of radius Ri located at n with the translational 
velocity Vi and the angular velocity This grain in- 
teracts with another grain j whenever overlapped; i.e. 
|ry| < Ri + Rj, where r,j = — rj. The interac- 
tion consists of two kinds of forces, each of which is 
normal and transverse to r^, respectively. Introduc- 
ing the unit normal vector ny = ry/|ry|, the normal 

(n) 

force acting on i, which is denoted by FJ - , is given by 
[/( £ y)+Cn.y ■ r y ] n y -, where e»j = l-\r i j\/(R i +R j ), A 
function /(e) describes elastic repulsion between grains. 
Here we test two models: /(e) = fee (the linear force) 
and /(e) = fee 3 / 2 (the Hertzian force) Note that 

the constant k/d 2 is on the order of the Young's mod- 
ulus of the grains. In order to define the transverse 
force, we utilize the relative tangential velocity v.y de- 
fined by (tij - riij • iij) + (RiCti + Rjftj)/(Ri + Rj) x r y 
and introduce the relative tangential displacement vec- 
tor A-^ = J roll divj* . The subscript in the integral 
indicates that the integral is performed when the con- 
tact is rolling; i.e., |A^| < k t or Ay ■ vfj < 0. Then 
the tangential force acting on the particle i is written 
as — min(/xrjy/|rjj|, fc t Ay )|Fy |. In the case that (J, = 0, 
the tangential force vanishes and the rotation of particles 
does not affect the translational motion. The parameter 
values adopted in the present simulation are given in TA- 
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TABLE I: The parameters of the discrete element simulation. 



polydispersity i^^Jd/km k t d fi Pd 2 jk 

30 % 1 0.005 - 0.6 3.8 x 10" 5 - 1.1 x 10" 2 



uniform shear is realized. In this case, the shear rate 
is proportional to the sliding velocity of the walls; i.e., 
7 = V/L z . 

We investigate the behaviors of the friction coefficient 
of the system, F y /PS = M. The control parameters 
that affect the friction coefficient are the shear rate 7 
and the pressure P. It is useful to represent the control 
parameters in terms of nondimensional numbers, because 
the friction coefficient is a nondimensional number and 
hence must be a function of nondimensional numbers. 
Thus 7 and P are recast in the following forms; I = 
jy/m/Pd and II = Pd 2 /k. In particular, the former is 
referred to as the inertial number [16j . which dominates 
the frictional behavior of granular flows. Hereafter we 
discuss the nature of friction taking advantage of these 
nondimensional numbers. 

In order to grasp the main point of our result, it is 
convenient to begin with the frictionless particles; i.e., 
fi = 0. Recall that we test two models, each of which has 
different interaction: the Hertzian contact model and the 
linear force model. As shown in FIG. [U friction coeffi- 
cients of these two models are collapsed on the following 
master curve. 

M = M + si*, (1) 

where Mo denotes the friction coefficient for 7^0. Here 
Mq ~ 0.06. Note that the effect of the inertial number 
is expressed by a power law, si* , where cj> = 0.27 ± 0.05. 
The prefactor s is approximately 0.37 in the linear force 
model, while it is somewhat larger (s ~ 0.45) in the 
Hertzian contact model. 

In order to check the universality of Eq. , we wish 
to confirm independence of our results on the details of 
the model. First we discuss the effect of the tangential 
force between particles. In FIG. shown are the friction 
coefficients of the models in which fi = 0.2 and 0.6. It is 
noteworthy that they range from 0.3 to 0.4, which are not 
significantly discrepant from those obtained by an exper- 
iment on spherical glass beads [l8|. More importantly, 
the friction coefficients again obey Eq. ([I]) with <f> — 0.3 
regardless of the value of \i and the force model (the lin- 
ear or the Hertzian). Indeed the friction coefficients of 
the both models are almost the same. We also remark 
that the factor s does not depend on /1; s = 0.33 ± 0.03 
for fi = 0, 0.2, and 0.6. 

On the other hand, Mq depends on /1. In the linear 
force model, M ~ 0.06 for fi = 0, while M a ~ 0.26 for 
/i = 0.2, and M ~ 0.4 for u = 0.6. Similar dependence 
was also observed in Refs. [HI, [111 . Note that Mq also de- 
pends on the force model in the case that \i = 0; Mq ~ 
in the Hertzian model while Mo is not negligible in the 
linear model. Although Mq looks like the static friction 
coefficient, note that it is defined in the 7^0 limit and 
different from the static friction coefficient, above which 
a static system begins to flow. In order to distinguish the 
two concepts, Mo is referred to as dynamic yield strength. 
The difference is important when we consider the stabil- 
ity of slip, as will be discussed in the last paragraph of 
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The configuration of the system mimics a typical ex- 
periment on granular layers subjected to simple shear. 
Note that there is no gravity in the system. The system 
spans L x x L y x L z volume, and is periodic in the x and 
the y directions. We prepare two systems of different as- 
pect ratio, each of which contains approximately 10,000 
particles: 25c? x 25g! x 80!, and 15d x 15d x 25d. As we 
shall discussed later, the difference of the aspect ratio 
does not affect the rheology. In the z direction, there 
exist two rough walls that consist of the same kind of 
particles as those in the bulk. The particles that consist 
walls are randomly placed on the boundary and their rel- 
ative positions are fixed. The walls are parallel to each 
other and displaced antiparallel along the y axis at con- 
stant velocities ±V/2, while they are prohibited to move 
along the x axis. One of the walls is allowed to move 
along the z axis so that the pressure is kept constant at 
P. Using the mass of the wall M w that is defined as 
the sum of the masses of the constituent particles, the z 
coordinate of the wall Z w is described by the following 
equation of motion; M W Z W — F z — PS, where F denotes 
the sum of the forces between the wall particles and the 
bulk particles, and S denotes the area of the wall. Then 
the z component of the velocity of the wall particles is 
given by Z w . Note that the friction coefficient of the 
system is defined by F y /PS. 

The system reaches a steady state after a certain 
amount of displacement of the walls. We judge that the 
system reaches a steady state if each of the following 
quantities does not show apparent trends and seems to 
fluctuate around a certain value. The monitored quan- 
tities are the friction coefficient, the z coordinate of the 
wall (i.e., the density), and the granular temperature. 
Also snapshots of the velocity profile are observed to en- 
sure the realization of uniform shear flow. We confirm 
that the transient behaviors of the friction coefficient and 
of the volume increase are quite similar to those observed 
in experiments. Here we do not investigate such tran- 
sients and restrict ourselves to steady-state friction. 

Because uniform shear flow is unstable in a certain 
class of granular systems, we must check the internal ve- 
locity profiles at steady states. There is a strict tendency 
that shear flow is localized near the walls in the case that 
the confining pressure is small and/or the sliding veloc- 
ity is large. This kind of spatial inhomogencity is rather 
ubiquitous in granular flow, and is extensively investi- 
gated [bH . Il5| . In our simulation, uniform shear flow is 
realized at lower sliding velocities and higher confining 
pressures. Here we discuss exclusively the case in which 
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FIG. 1: The friction coefficients M in the models without 
tangential force, i.e., if — 0. The horizontal axis / denotes 
the nondimensional shear rates, V / L^Jm/Pd. The shape of 
the symbols and the confining pressure are in one-to-one cor- 
respondence: the squares to IT = 3.8 x 10 -5 , the circles to 
II = 1.9 x 10" 3 , the triangles to II = 3.9 x 1CT 3 , and the 
diamonds to II = 1.1 x 1CP 2 . The layer thickness L z /d ~ 8 
for the squares and the circles, while L z /d ~ 26 for the trian- 
gles and the diamonds. The solid lines denote Eq. fTJ with 
4> — 0.3. (a) Friction coefficient of the Hertzian force model, 
(b) Friction coefficient of the linear force model. 



this Rapid Communication. 

It is important to notice that the data at different con- 
fining pressures collapse on the same curve by virtue of 
the inertial number. This suggests that the friction co- 
efficient of a dense granular material does not depend 
on IT, as long as it is small (in the present simulation 
3.8 x 1CT 5 < n < 1.1 x 1(T 2 ). Indeed, da Cruz et al. 
pd| found that the friction coefficient is independent of 
II in their notation) up to II < 2.5 x 1CP 2 in a two 
dimensional system. Therefore the independence on II 
is very likely within the accuracy of these simulations. 
While one can still expect that the dependence may ap- 
pear for larger II, such a case that II ~ 0.1 is meaningless 
as a model of a granular material. Note that II roughly 
corresponds to the average overlap length of contacts 
divided by the particle diameter; namely, the average 
strain of individual particles. 

Then we discuss the effect of inelasticity that is mod- 
eled by the viscous coefficient (. The co rrespo nding 
nondimensional number £ is defined by Qydjkm. We 
find that decrease of £ reduces the friction coefficient in 
the region where I > 0.01, while the frictional strength 
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FIG. 2: The friction coefficients in the models with tan- 
gential force (/i = 0.2, 0.6). The shape of the symbols and 
the confining pressure are in one-to-one correspondence as in 
FIG. Q] The blank symbols denote the friction coefficient of 
the linear force model with if = 0.2, while the symbols of 
vertically-striped pattern denote that of /i = 0.6. The cir- 
cles of horizontally-striped pattern denote the Hertzian force 
model with [i = 0.2. The lines denote Eq. (JTJ with </> = 0.3. 

is independent of £ for smaller / region. This behavior is 
consistent with those obtained in Refs. [HI, HH. Never- 
theless, it can be still described by Eq. ([TJ) with s being 
a smaller value. For example, the friction coefficient of 
a system in which £ = 0.05 is described by s ~ 0.27 
with almost the same values of M and S. However, the 
functional form of s(£) is not clear at this point. 

From the discussions so far, we can conclude that the 
details of the present model do not affect the validity 
of Eq. ([T]), which is the main result of this study. Im- 
portantly, the exponent <f> seems to be universal; it is ap- 
proximately 0.3 regardless of the details of the model and 
the control parameters. The velocity-strengthening na- 
ture of this friction law does not contradict experiments 
0) S @, S Hi • br addition, it illustrates universality of 
power-law in rheological properties of random media [2(| 
including foams plf and human neutrophils [12]. In the 
following we discuss four important points that are pe- 
ripherally related to the main result. 

First, we discuss the dependence of the volume frac- 
tion to the inertial number. Surprisingly, decrease of the 
volume fraction caused by shear flow is also described by 
a power-law. 

V Q - V = 821 s , (2) 

where vq is the volume fraction in the 7 — > limit. Note 
that the constants S2 and 6 do not depend on the details 
of the model. Figure [3] shows that all of the data obtained 
in our model collapse on Eq. © with S2 ~ 0.11 and 5 = 
0.56 ± 0.02. This dilatation law also illustrates ubiquity 
of power-law in granular materials. 

The next point we wish to discuss is the relation be- 
tween the present result and power-law rhcology in sys- 
tems at constant volume. In particular, Xu and O'Hern 
[23| found a power-law relation between the shear stress 
and the shear rate in a two dimensional granular mate- 
rial consisting of frictionless particles. They estimated 
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FIG. 3: The dilatation law. Decrease of the volume fraction 
Av — i/o — v is plotted as a function of the inertial number. 
Note that vo is the volume fraction in the 7 — » limit, which 
is estimated by extrapolation. The symbol legends are the 
same as those in FIGS. Q] and The line denotes Eq. 

the exponent to be 0.65. However, at constant volume 
condition, the pressure also depends on the shear rate so 
that the behavior of the friction coefficient is generally 
different from that of the shear stress. Therefore, power- 
law friction in systems under constant volume condition 
are not directly related to the present result. See Ref. 
[24j for more detailed discussions on this subject. 

The third point we wish to discuss is the effect of di- 
mensionality. In contrast to the present study, da Cruz 
et al. [ll.] obtained a linear friction law in a two dimen- 
sional system. The difference may be attributed to the 
dimensionality of the systems, which affect the nature 



of contacts between particles. In particular, the angular 
distribution of the tangential force is strongly anisotropic 
in two dimensional systems, while such anisotropy is not 
observed in our three dimensional system. Accordingly, 
in the case of frictionless particles, their system exhibited 
a friction law that is quite similar to ours. 

As the fourth point of interest, we discuss relevance of 
our result to earthquake mechanics by comparing it to 
a friction law recently proposed by Jop et al. 0], which 
seems to be validated in experiments on inclined plane 
flow. We stress that such flow is characterized by rela- 
tively large inertial number (typically / > 10 -1 ), while 
our simulation involves much smaller / values (/ > 10~ 4 ) 
as shown in FIGS. [T] and O In short, Eq. involves 
much smaller I region than Jop et al. have investigated. 

Such small inertial numbers correspond to a typical 
configuration of seismic motion of faults. For example, 
in the case that d = 1 mm, V = 1 m/s, L z = 4 cm, 
and P = 100 MPa, the corresponding inertial number is 
10~ 4 . However, one may wonder that the friction law 
Eq. (TTJ) cannot lead to stick-slip motion of faults because 
the friction law found here is velocity-strengthening. Re- 
call that we discuss exclusively stationary-state dynamic 
friction. Taking static friction into account, unstable 
slip is inevitable because static friction is always stronger 
than dynamic friction, which is mainly due to dilatation. 
Therefore power-law friction in stationary states does not 
contradict unstable slip on faults. 
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